#include "findMisMatch_gpu.h"

__global__ void baseMisMatch(char *d_read, char *d_ref, unsigned int *d_readPos, unsigned int *d_refPos, unsigned char* d_numMisMatch, unsigned char** d_misPos, unsigned char** d_misBase, const int readLength, unsigned int numSA)
{
        __shared__ unsigned int misMatch[FQUERY_PER_BLOCK];
	unsigned int x = threadIdx.x;
	unsigned int bOffset = blockIdx.x*FQUERY_PER_BLOCK  + x;
	unsigned int refPos;
	unsigned int readPos;
        unsigned int tOffset;
	while(bOffset < numSA)
	{
		__syncthreads();
		refPos = d_refPos[bOffset];
		readPos = d_readPos[bOffset];
                __syncthreads();
		misMatch[x] = 0;
		for(tOffset = threadIdx.y; tOffset < readLength; tOffset += FWARP_SIZE)
	        {        
			if(d_ref[refPos+tOffset] != d_read[readPos+tOffset])
	                {      
   			        unsigned int old =  (atomicInc(&misMatch[x], FWARP_SIZE))%2;
		                d_misPos[old][bOffset] = tOffset;
		                d_misBase[old][bOffset] = d_read[readPos+tOffset];
		     
	                }
               }
	       __syncthreads();
	       d_numMisMatch[bOffset] = misMatch[x];
               bOffset += FTHREAD_NUM;
	}
}

__global__ void baseMisMatch_v1(char *d_read, char *d_ref, unsigned int *d_readPos, unsigned int *d_refPos, unsigned char* d_numMisMatch, unsigned char** d_misPos, unsigned char** d_misBase, const int readLength, unsigned int numSA, unsigned offset)
{
        __shared__ unsigned int misMatch[FQUERY_PER_BLOCK];
	unsigned int x = threadIdx.x;
	unsigned int bOffset = blockIdx.x*FQUERY_PER_BLOCK  + x + offset;
	if(bOffset >= numSA)
		return;
        unsigned int tOffset = threadIdx.y;
	unsigned int refPos = d_refPos[bOffset];
	unsigned int readPos = d_readPos[bOffset];
	__syncthreads();
	misMatch[x] = 0;
	while(tOffset < readLength)
        {        
		if(d_ref[refPos+tOffset] != d_read[readPos+tOffset])
	       {      
   		        unsigned int old =  (atomicInc(&misMatch[x], FWARP_SIZE))%2;
	                d_misPos[old][bOffset] = tOffset;
	                d_misBase[old][bOffset] = d_read[readPos+tOffset];
	     
	       }
	       tOffset += FWARP_SIZE;
        }
	__syncthreads();
        d_numMisMatch[bOffset] = misMatch[x];
}
__global__ void multiplyLength(unsigned int* d_readPos, unsigned int readLength, unsigned int numSA)
{
	unsigned int offset = blockIdx.x*FBLOCK_SIZE + threadIdx.x;
	while(offset < numSA)
	{
		d_readPos[offset] *= readLength;
		offset += FCALPOS_THREAD_NUM;
	}
}


void handleError(cudaError_t err, int line)
{
	if (err != cudaSuccess)
		printf("line number: %d  %s\n", line, cudaGetErrorString(err));
}
void calPos(unsigned int*offset, unsigned int *pos, unsigned int *d_readPos, unsigned int *d_refPos, unsigned int numSA, unsigned int readLength)
{
	unsigned int i;
        unsigned int* tempPos;
	cudaHostAlloc(&tempPos, UINT_SIZE*numSA, cudaHostAllocDefault);
	for(i = 0; i < numSA; i++)
		tempPos[i] = pos[offset[i]];
	HANDLE_ERROR(cudaMemcpy(d_refPos, tempPos, UINT_SIZE*numSA, cudaMemcpyHostToDevice));
	multiplyLength<<<FBLOCK_NUM, FBLOCK_SIZE>>>(d_readPos, readLength, numSA);
	
	cudaThreadSynchronize();
	HANDLE_ERROR(cudaFreeHost(tempPos));
}

void findMisMatch(char* read, char *ref, unsigned int *queryId, unsigned int *pos, unsigned int *offset, Info *info)        
{	
	char *d_read;
	char *d_ref;
	unsigned char **d_misPos;
	unsigned char **d_misBase;
	unsigned char *tempBase[2];
	unsigned char *tempPos[2];
	unsigned char *misBase1 = info->misBase1;
	unsigned char *misBase2 = info->misBase2;
	unsigned char *misPos1 = info->misPos1;
	unsigned char *misPos2 = info->misPos2;
	unsigned char *numMisMatch = info->numMisMatch;
	unsigned int *d_readPos;
	unsigned int *d_refPos;
	unsigned int readLen = info->readStringLength;
	unsigned int refLen = info->refStringLength;
	unsigned int readLength = info->readLength;
	unsigned int numSA = info->numSA;
	unsigned char* d_numMisMatch;
	unsigned char* d_misBase1;
	unsigned char* d_misBase2;
	unsigned char* d_misPos1;
	unsigned char* d_misPos2;
	
	size_t biCharPoints = sizeof(unsigned char*)*2;
	dim3 blockSize(FQUERY_PER_BLOCK, FWARP_SIZE);
        float calTime = 0;
	float sumTime = 0;
	float memTime = 0;
	float calPTime = 0;

	timeval memStart, memEnd;
	timeval calStart, calEnd;
	timeval sumStart, sumEnd;
	timeval calPStart, calPEnd;
	gettimeofday(&memStart, NULL);
	gettimeofday(&sumStart, NULL);

	printf("findMisMatch: processing...\n");
	cudaFuncSetCacheConfig(baseMisMatch, cudaFuncCachePreferL1);
	cudaFuncSetCacheConfig(multiplyLength, cudaFuncCachePreferL1);
	
	HANDLE_ERROR(cudaMalloc(&d_read, readLen));
	HANDLE_ERROR(cudaMalloc(&d_ref, refLen));
	HANDLE_ERROR(cudaMalloc(&d_readPos, numSA*UINT_SIZE));
        HANDLE_ERROR(cudaMalloc(&d_refPos, numSA*UINT_SIZE));
	
	HANDLE_ERROR(cudaMalloc(&d_numMisMatch, numSA));
	HANDLE_ERROR(cudaMalloc(&d_misBase1, numSA));
	HANDLE_ERROR(cudaMalloc(&d_misBase2, numSA));
	HANDLE_ERROR(cudaMalloc(&d_misPos1, numSA));
	HANDLE_ERROR(cudaMalloc(&d_misPos2, numSA));
	
	HANDLE_ERROR(cudaMemset(d_numMisMatch,0, numSA));
	HANDLE_ERROR(cudaMemset(d_misBase1,0, numSA));
	HANDLE_ERROR(cudaMemset(d_misBase2,0, numSA));
	HANDLE_ERROR(cudaMemset(d_misPos1,0, numSA));
	HANDLE_ERROR(cudaMemset(d_misPos2,0, numSA));

	HANDLE_ERROR(cudaMalloc(&d_misPos,biCharPoints));
	HANDLE_ERROR(cudaMalloc(&d_misBase,biCharPoints));

	tempBase[0] = d_misBase1;
	tempBase[1] = d_misBase2;
	tempPos[0] = d_misPos1;
	tempPos[1] = d_misPos2;

	HANDLE_ERROR(cudaMemcpy(d_misPos, &tempPos, biCharPoints, cudaMemcpyHostToDevice));
	HANDLE_ERROR(cudaMemcpy(d_misBase, &tempBase, biCharPoints, cudaMemcpyHostToDevice));
	HANDLE_ERROR(cudaMemcpy(d_read, read, readLen, cudaMemcpyHostToDevice));
	HANDLE_ERROR(cudaMemcpy(d_ref, ref, refLen, cudaMemcpyHostToDevice));
	HANDLE_ERROR(cudaMemcpy(d_readPos, queryId, numSA*UINT_SIZE, cudaMemcpyHostToDevice));

	cudaThreadSynchronize();
	
	gettimeofday(&memEnd, NULL);
        memTime = (float)(FCLK*(memEnd.tv_sec - memStart.tv_sec)+(memEnd.tv_usec - memStart.tv_usec))/FCLK;
	
	gettimeofday(&calPStart, NULL);
	calPos(offset, pos, d_readPos, d_refPos, numSA, readLength);
	gettimeofday(&calPEnd, NULL);
	
	gettimeofday(&calStart, NULL);
	for(unsigned i = 0; i <= numSA/FTHREAD_NUM; i++)
	{
		baseMisMatch_v1<<<FBLOCK_NUM, blockSize>>>(d_read, d_ref, d_readPos, d_refPos, d_numMisMatch, d_misPos, d_misBase, readLength, numSA, i*FTHREAD_NUM);	
		cutilCheckMsg("baseMisMatch_v1");
	}
	cutilSafeCall(cudaThreadSynchronize()); 
	gettimeofday(&calEnd, NULL);
	gettimeofday(&memStart, NULL);

	HANDLE_ERROR(cudaMemcpy(misBase1, d_misBase1, numSA, cudaMemcpyDeviceToHost));
	HANDLE_ERROR(cudaMemcpy(misBase2, d_misBase2, numSA, cudaMemcpyDeviceToHost));
	HANDLE_ERROR(cudaMemcpy(misPos1, d_misPos1, numSA, cudaMemcpyDeviceToHost));
	HANDLE_ERROR(cudaMemcpy(misPos2, d_misPos2, numSA, cudaMemcpyDeviceToHost));
	HANDLE_ERROR(cudaMemcpy(numMisMatch, d_numMisMatch, numSA, cudaMemcpyDeviceToHost));
	
	cutilSafeCall(cudaThreadSynchronize()); 
/*	
	for(int i=0; i<numSA; i++)
		printf("%c ",misBase1[i]);
	printf("\n\n");
	for(int i=0; i<numSA; i++)
		printf("%d ",misPos1[i]);
	printf("\n\n");
	for(int i=0; i<numSA; i++)
		printf("%c ",misBase2[i]);
	printf("\n\n");
	for(int i=0; i<numSA; i++)
		printf("%d ",misPos2[i]);
	printf("\n\n");
*/
	
	HANDLE_ERROR(cudaFree(d_read));
	HANDLE_ERROR(cudaFree(d_ref));
	HANDLE_ERROR(cudaFree(d_readPos));
	HANDLE_ERROR(cudaFree(d_refPos));
	HANDLE_ERROR(cudaFree(d_numMisMatch));
	HANDLE_ERROR(cudaFree(d_misBase1));
	HANDLE_ERROR(cudaFree(d_misBase2));
	HANDLE_ERROR(cudaFree(d_misPos1));
	HANDLE_ERROR(cudaFree(d_misPos2));

	cutilSafeCall(cudaThreadSynchronize()); 
	
	gettimeofday(&memEnd, NULL);
	gettimeofday(&sumEnd, NULL);

	calTime = (float)(FCLK*(calEnd.tv_sec - calStart.tv_sec)+(calEnd.tv_usec - calStart.tv_usec))/FCLK;
	sumTime = (float)(FCLK*(sumEnd.tv_sec - sumStart.tv_sec)+(sumEnd.tv_usec - sumStart.tv_usec))/FCLK;
	memTime += (float)(FCLK*(memEnd.tv_sec - memStart.tv_sec)+(memEnd.tv_usec - memStart.tv_usec))/FCLK;
	calPTime = (float)(FCLK*(calPEnd.tv_sec - calPStart.tv_sec)+(calPEnd.tv_usec - calPStart.tv_usec))/FCLK;
	
	printf("findMisMatch: finished\n");
	printf("findMisMatch: Total time: %.2f s\n", sumTime);
	printf("findMisMatch: Time used copying data between Host and Device: %.2f s\n", memTime);
	printf("findMisMatch: Time used calculating readPos and refPos: %.2f s\n", calPTime);
	printf("findMisMatch: Calculation time: %.2f s\n", calTime);
}
